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ABSTRACT 


The transient phase change process that occurs when a vapor contacts a very cold 
inclined plate is analyzed. The vapor condenses to form a growing liquid layer that 
flows down the plate. Below this two-phase region there is a three-phase region where 
the thickening condensate layer flows over a growing frozen layer. Heat and mass are 
transferred across moving phase boundaries. The solution involves the Karmen- 
Pohlhausen method which results in quasi -linear partial differential equations for the 
phase thicknesses that are solved by characteristic and finite difference methods. Nu- 
merous two-phase special cases are also discussed. 
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THREE-PHASE HEAT TRANSFER: TRANSIENT CONDENSING WITH 
FREEZING ON A NONISOTHERMAL INCLINED PLATE 
by William A. Olsen, Jr., and Frank B. Molls 
Lewis Research Center 

SUMMARY 

The subject of this analytical study is transient condensing with freezing of a slowly 
moving saturated vapor upon an inclined plate. The plate has no thermal capacity and is 
cooled by a fluid that is below the freezing point of the vapor. At the top of the plate, 
there is a two-phase region where a growing condensate liquid layer flows (laminar 
flow) down the plate. Below that region, there is a three-phase region where the grow- 
ing condensate layer flows over a thickening frozen layer. 

The solution to the problem employs the Karman-Pohlhausen integral method with 
assumed linear temperature profiles across the phase layers. This method reduces the 
problem to quasi-linear first-order hyperbolic partial differential equations, with the 
phase thicknesses as dependent variables. For the two-phase region, the solution in- 
volves a single characteristic curve, which is represented by algebraic equations. 
Families of characteristic curves, described by a network of ordinary differential equa- 
tions, are involved in the solution for the three-phase region. The hyperbolic partial 
differential equations for the three-phase region are also solved by a direct finite dif- 
ference equation approximation for comparison purposes . 

Results of the analysis indicate that freezing will start at and beyond a fixed point, 
as soon as the insulating condensate layer is thick enough for the wall temperature at 
that point to drop to the freezing temperature. The solid and condensate layers, which 
have a shape similar to boundary layers, grow until they attain a steady-state profile. 

This analysis specializes to a number of cases which are discussed herein. Some 
of these cases are not covered in the literature. 



INTRODUCTION 


This report deals with the transient and steady- state phase change processes that 
occur when a vapor suddenly comes in contact with an inclined nonisothermal plate. The 
plate is cooled by a fluid that is below the freezing point of the vapor. In one area of 
the plate condensing occurs, while in another area condensing accompanied by freezing 
occurs. The literature for two-phase problems (e.g. , freezing, melting, or condensing) 
is extensive and is summarized in references 1 and 2. However, the three-phase heat- 
transfer problem, which combines two-phase problems, has apparently received little 
attention. There has been little interest in the three-phase problem because it occurs 
only rarely in nature, and because commercial applications have been lacking. The 
reason for this is that large temperature differences are normally required in order to 
span across the thermodynamic phase space from the vapor to the solid region by cross- 
ing the liquid region. Now however, large temperature differences can often occur be- 
cause of the cold environments of space missions, and because cryogenic liquids and 
liquid metals are broadly applied. A major problem of heat-transfer systems that oper- 
ate in an environment colder than the freezing point of the vapor is that a solid layer 
could occur in the vapor flow passages. Any such solid layer that forms there would 
usually be detrimental in that it restricts flow or reduces the heat transfer. 

In reference 3 the transient horizontal plate and steady-state inclined plate cases of 
the three-phase problem were considered. The problem considered herein in an exten- 
sion to the work of reference 3. This extension considers transient condensing and 
freezing of a slowly moving vapor upon a nonisothermal inclined plate. The condensate 
is subjected to a constant body force (e.g. , an inclined plate in a constant gravity field). 
The plate shown in figure 1 is cooled by a coolant of constant temperature and constant 
heat-transfer coefficient (e. g. , a well-mixed coolant bath). Initially the plate is dry, 
and suddenly conditions change such that condensing from the vapor commences. This 
may later be followed by freezing of the condensate liquid beyond some point on the plate 
(X on fig. 1). Before steady- state phase thicknesses are attained, there is a growing 
condensate layer (two-phase region) above X„. Below X there is a three-phase re- 
gion, where the condensate layer flows over a solid layer that grows by freezing some of 
the condensate liquid. Heat and mass are transferred across the moving phase bound- 
aries in both the two- and three-phase regions. 

In order to continue the description of the physical model, it is necessary to men- 
tion two effects of a noncondensable gas. The experiment of reference 3 dealt with con- 
densing and freezing of water vapor, in the presence of some air, upon a liquid- 
nitrogen- cooled vertical cold plate. These experimental results qualitatively indicated 
that a noncondensable gas (air) has a large effect upon the manner of formation and the 
type of solid layer formed from the vapor (water vapor). As the noncondensable gas 
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fraction increased, the type of solid formed from the vapor changed from an ice-like 
solid to an ice-and-frost-laminated composite solid, and finally to frost alone at large 
air fractions. Only the first case, condensing from a pure vapor where a solid (ice-like) 
layer is formed, is considered in the transient case of this report. The second effect is 
that a small amount of non condensable gas (small enough for an ice -like solid) offers 
considerable thermal resistance to heat transfer such that a much large solid layer 
grows than with a pure vapor. This second effect is quantitatively considered in appen- 
dix C by a modification to the limiting case of steady state. 

An exact solution (i. e. , algebraic equations) to this problem is unlikely. The ap- 
proximate Karman-Pohlhausen integral method was recommended for three-phase prob- 
lems in reference 3, where many methods were compared. Accordingly, this is used 
herein and results in a set of quasi-linear first-order partial differential equations to 
solve. Both characteristics and finite difference techniques are employed in their solu- 
tion. The accuracy of these methods is compared to the exact solutions of special 
cases. 

The analysis of transient condensing with freezing specializes to a number of cases, 
some of which have not been reported in the literature. Transient condensing and freez- 
ing on a horizontal plate and steady- state condensing and freezing on an inclined noniso- 
thermal plate are two such cases. Transient condensing on a nonisothermal plate is an- 
other. With minor changes in the initial conditions, the analysis can handle transient 
condensing with melting of a thin slab of solid by its vapor; and, with a minor modifica- 
tion, transient freezing of a liquid flowing past a cold plate can be analyzed. 


DERIVATION OF DIFFERENTIAL EQUATIONS 

Consider the plate of no thermal capacity (i. e. , a thin metal plate a short time after 
the start of the transient) shown in figure 1. One side of the plate is in contact with a 
coolant of constant heat-transfer coefficient h c and constant temperature T . Sud- 
denly, a slowly moving vapor comes in uniform contact with the opposite side of the cold 
plate and the vapor starts to condense while the pressure P is held constant (fig. 1(a)). 
The condensate liquid layer will grow thicker while it flows down the plate (fig. 1(b)). 

The flow is caused by some steady body force (e. g. , gravity with the plate inclined from 
the horizontal) and/or vapor drag caused by momentum transfer and friction. As the 
insulating layer of condensate liquid thickens, the temperature of the plate surface in 
contact with the phases T w will fall. The coolant temperature T c may be below the 
freezing temperature of the vapor Tg^. If it is, may fall to the freezing tempera- 
ture T SL* When this happens, a solidified layer will start to grow at those areas of the 
plate. Further growth of the condensate liquid will result in a corresponding growth in 
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the solidified layer until a steady state is attained (fig. 1(c)). Because of the finite ther- 
mal resistance of the coolant and plate b, the start of freezing will occur some distance, 
x = X g (t) , below the start of the cold plate, x = 0. This problem, therefore, consists 

of two regions: one where transient condensing (two phases) occurs for x < X , and an- 

s 

other where transient condensing and freezing (three phases) occur for x > X . Should 

s 

the coolant and plate thermal resistance b be zero, the plate would be at a constant 
wall temperature, T = T , and freezing would start at x = 0. 

W L 

This three-phase problem with its moving phase boundaries can be quite compli- 
cated, such that a number of simplifying assumptions that do not greatly compromise 
the physics of the problem are necessary. The three-phase problem is essentially a 
combination of the two-phase problems of condensing and freezing. The literature for 
each of these two-phase studies is extensive. These two-phase studies are used herein 
to supply physically reasonable simplifying assumptions for the three-phase problem. 

Consider the condensate liquid layer first. The literature for steady-state condens- 
ing on a vertical isothermal plate is well summarized in references 2 and 4, and trans- 
ent condensing on a vertical isothermal plate has been studied in reference 5. These 
references indicate that the physical assumptions necessary to obtain the equations of 
Nusselt's original simple model for condensing on a plate, as described in reference 6, 
are adequate for condensing of a pure slowly moving vapor. Furthermore, the Nusselt 
model can be adapted to successfully describe the effect of noncondensable gas and also 
liquid-metal condensing. Based upon the simple Nusselt model, the following assump- 
tions are made for the analysis of the condensate layer in this report: no waves, small 
phase and wall boundary curvature, constant properties, axial gradients are much less 
than normal gradients and axial velocities are much larger than normal velocities 
(boundary layer assumptions), saturated pure vapor at constant pressure, laminar flow, 
and, finally, the momentum convection terms (e. g. , 3u^/3t + u^ 3u^/9x + v^ 3u^/3y) have 
a negligible effect upon the flow in the thin condensate layer. Unlike the Nusselt model, 
the condensate energy convection terms are included herein; and the wall temperature 
is allowed to vary, but axial heat transfer is neglected. The coolant temperature and 
the thermal resistance of the coolant and plate are assumed to be constant. 

By these assumptions, the general continuity, momentum, and energy equations that 
are summarized in reference 7 are sim plied to the following equations that describe the 
condensate liquid layer. The coordinate system is fixed to the plate, as shown in fig- 
ure 1. Continuity is described by 


d A 

3t 


+ — (P,u 7 )+— (p,v ) = 0 

3x L L 3y 1 c 


( 1 ) 
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where the constant density assumption has not as yet been applied to equation (1) for con- 
venience. The momentum equation for the thin slowly moving condensate layer simpli- 
fies greatly and involves only a bouyant force term and a viscous drag term. 


(Pi ~ P y )S x + 



( 2 ) 


The energy equation for the condensate layer is given by equation (3). 


Pl c l 



+ «7 


dx 


+ V 7 



8 ^ 

3y 



(3) 


Equations (1) to (3) are also used in reference 5 to describe transient condensing on an 
isothermal vertical plate. 

Just as with the condensate layer, it can be expected that the axial gradients in a thin 
stationary solid layer are small compared to the normal gradients. By making the addi- 
tional assumptions that heat is transfered in the solid solely by conduction and that the 
solid properties are constant, the following equation can be written to describe heat 
transfer in that layer: 


9T 


Ps c 


s s 


at 


% k 

3 3y 2 


(4) 


Reference 1 used this equation for one-dimensional transient freezing on a plate. 

The boundary conditions of the variables of the problem (e. g. , T(x, y, t) , 

Y L v( x >t), q(x, y, t) , u(x, y, t), etc. ) are listed belcw. The term T(x, y, t) means 
that the temperature is a function of x, y, and t. 

(1) Initially, the slowly moving saturated vapor is in uniform contact with a cold dry 
plate so that for all x at t = 0 the following are the boundary conditions: 


^LV^ X ’ = 0 

(5) 

^Sl/ X ’ = o 

(6) 

T w< x ’ °> = T c 

0 ) 
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(2) At the leading edge of the cold plate (x = 0), the condensate thickness is zero. 

Y LV < 0 ,t>= 0 (8) 

(3) Along the inner liquid condensate layer boundary (fig. 1), there are two regions 
where the boundary conditions differ. 

(a) Along the plate (y = 0) where there is no solid layer (x < X (t)), the fol- 

s 

lowing conditions occur: 


Uj(x,y=0,t) = 0 
Tj(x,0,t) =T w (x,t) 


(9) 

( 10 ) 


-q z (x, 0,t> 




(x, 0, t) 



( 11 ) 


where b is the combined coolant and plate thermal resistance, b = 1/h + d m /k m . 

(b) At the point where freezing starts on the plate (x = X (t)), the wall tern- 

b 

perature is at the freezing point. 


Y SL( x = X s< t >’ t > =0 

(12) 

T (X , t) = T ct 
w x s’ 7 SL 

(13) 


(c) The boundary conditions along y = Yg^(x, t) (the solid-liquid interface), 

which apply where there is a solid layer at x > X (t) , are given by 

s 

-q s (x,y=Y SL ,t> +q z (x,Y SL ,t) 


3T 


ST, 


- k g — -2. (x, Y gL , t) - k^ —± (x, Y gL , t) 


•(r) SL ^-\ A / fr ^ 


(14) 
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U Z^ X,Y SL’ 9 “ 0 


(15) 


T Z^ X> Y SL’^ - T s ( x ’ Y SL’ ^ - t sl 


(16), (17) 


(4) Along the plate (y = 0) where there is a solid layer (x > X g (t)), the following 
conditions are used: 


T s <x,y=0,t> = T w (x,t> 


(18) 


-q (x, 0, t) = k — (x, 0, t) = Tw ^ X ’ - — = /Q\ ( x , t) 
3y b \aI 


(19) 


(5) Along the entire liquid-vapor interface (x > 0) which follows y = Y Ly (x, t) , the 
boundary conditions are 


Vx, Y L V >t> - Tj j y(P) (20) 

"V x ’ Y LV’t> = k £ ~ (x,Y LV ,t) = -0^ L LV =^— ^ L lv (21) 

where the vapor was assumed to be saturated at constant pressure so that 
0T v /0y(x, Y lv , t) = 0. The basis for equations (13), (14), (16), (17), (20), and (21) is 
discussed in appendix D of reference 3 and in reference 8. The final boundary condition 
is the shear stress at y = Ygj^(x, t). Reference 2, where this boundary condition is dis- 
cussed in detail, indicates that this shear stress is composed of contributions from the 
momentum transfer of condensation and form drag. Unfortunately, even low vapor ve- 
locities can readily cause considerable waviness of the liquid- vapor interface, which ap- 
preciably affects heat transfer. This situation adds too much complication to be con- 
sidered herein; therefore, it is assumed that the vapor velocity is low enough (e. g. , 

U y < 10 ft/sec (3 m/sec)) that the effect of vapor drag can be neglected. The effect of 
the small ripple waves caused by gravity is also neglected. The resulting shear stress 
boundary condition is given by equation (22). 

0U 7 

-i<x,Y LV ,t> ho (22) 

0y J - ,v 

This analysis is therefore realistic for external condensing (e. g. , a plate condenser) 
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but not for internal condensing (e. g. , condensing inside a tube), where the vapor velocity 
is very high. 

The next task is to apply the Karman-Pohlhausen integral method to the differential 
equations (eq. (1) to (4)) by first integrating them over the time- and position-varying 
phase thicknesses for a dx element. Following this operation, temperature profiles are 
assumed and the boundary conditions applied. These operations will result in a simpler 
system of partial differential equations, where the phase thicknesses are the dependent 
variables. 

A mass balance on the liquid condensate layer is obtained from the continuity equa- 
tion (eq. (1)) by this method (see appendix D in ref. 3 for more detail). Wherever there 

is a solid layer (x > X (t)), these operations result in the following integro-differential 

s 

equation. 



i 


Y LV < x > t > 


U 7 


'Y sl (x, t) 


dy +Pl 





where the constant density assumption has finally been applied. Wherever there is no 
solid layer (x < X g (t)), continuity of mass becomes 



f 

J 0 


LV 


<x,t> 


u i d y + Pi 


dY 


LV 


9t 



(23b) 


The mass balance on the solid layer is similarly derived. 


Ps 




(24) 


The energy equation for the condensate layer (eq. (3)) is similarly integrated. Where- 
ever there is a solid layer (x > X o (t)), equation (25) describes the energy balance. 
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pjC 


ri 


at 


pLV^*) 

J Y si/ x > 


(T z • T L y) d y + 


Pfy-O- r LV<X,t> a 

3 V Y SL< x . t > 


i z ( T z - T Ly )dy 


+ c z^ t lv - t sl) (f ) = - q z <*> Y LV > t) + V x > Y SL> *> 

' A/ SL 


( 25 ) 


Equation (26) describes the energy balance on the condensate layer where there is no 
solid layer at x< X g (t). 


p l c l~ 


at 


/ 


LV 


(x,t) 


(T * - T Lv> d y + pfi — 


jf 


LV 


<M> 


u ^(T^ - T j^y) dy q^(x, Y^y, t) 


+ q 2 (x, 0, t> (26) 

Integrating equation (4) over the solid layer thickness results in the following equation 
for the solid layer: 


p c„ 
's s 


a / ,Y sL< x ’ t > 


at 


X 


(T s - T 


SL 


)dy = 


q s (x,Y 


SL 


|5 t) + q (x,0,t) 


(27) 


The integration of the momentum equation (eq. (2)) is considered next. Because of 
the no-vapor-shear assumption, this equation can be very simply integrated. By using 
the velocity boundary conditions (eqs. (15) and (27)), the following velocity profile is ob- 
tained for where there is a solid layer at x > x <t>. 

S 


u i = V y > 


Y SL V 2 

2 “ Y + YLV ^ y ~ Ysl ^ 


And wherever there is no solid layer, x < X (t) 

s 



(28a) 


(28b) 


9 



where 


n 


(P; - P y )g x 


M 


l 


(28c) 


This velocity profile is now used in the mass and energy balances (eqs. (23), (25), 
and (26)) along with assumed temperature profiles across the phase layers in order to 
complete the indicated integrations. 

References 1 and 4 have shown that linear temperature profiles are good approxi- 
mations for the temperature distribution across thin phase layers where there is a phase 
change. With this physical insight, the following linear temperature profiles are as- 
sumed, which satisfy the temperature boundary conditions (e.g , eqs. (10), (16) to (18), 
and (20)). For the condensate layer where there is freezing (x > X (t)), assume the 
following linear temperature profile: 


T l - t lv * < t lv - m/ ' ^ - ') < 29 > 

\ LV ■ *SL / 

For the condensate layer where there is no freezing (x ^ X (t)), the following linear 

s 

relation is assumed: 


T i - t lv s < t lv - V x > (30) 

And for the solid layer at *>X s (t> 

T s - t sl = < t sl - *> (S [) 

Equations (30) and (31) involve a variable wall temperature, which is removed from 
the temperature profiles by using the boundary conditions. For x < X g (t), substitute 
equation (30) into equation (11) and solve for T (x, t) . 
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T w <x,t> = 


(32) 


bT i x 

Di LV + . A c 


LV 


+ b 


for x > X (t), substitute equation (31) into equation (19). 


bT + ^ x 
D1 SL + A c 


T w (x,t) = 


l SL 


+ b 


(33) 


These relations are used below to eliminate T w (x, t) when the integro -differential equa- 
tions are integrated and put into a form that contains only the phase thicknesses as vari- 
ables. 

The heat flux to the coolant (Q/Ai from the condensate layer at x< X (t) is de- 

C o 

termined by substituting equation (32) into equation (11). 



-q z (x,o,t> 



' T c 


+ b 


(34) 


The heat flux to the coolant from the solid layer at x > X g (t) is determined from equa- 
tions (33) and (19). 



-q s (x,o,t> 


T SL - T c 


l SL 


+ b 


(35) 


The heat flux into the solid layer from the condensate layer at x > X (t) is determined 

s 

from equations (29) and (14). 
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(36) 


-V^SL’*) = "V^SL^ 



¥ T LV - T SL> 

y lv “ y sl 



J SL 


And finally the heat flux into the liquid for all x is given by equation (21). 

The next steps in the analysis involve the evaluation of the integrals of equations (23) 
and (25) to (27) by using the velocity profiles, temperature profiles, and heat flux bound- 
ary conditions that have just been derived. These evaluations are greatly simplified by 
changing the Y L y and Yg L variables to the phase thickness variables 6 and A. 

Where there is no solid layer (x < X g (t)) 


y lv( x > 9 s 5 < x > t ) = 6 


(37) 


while for x>X g (t) 


Y L,v( x >t) " Y SL (x,t) - 6(x, t) - 6 


(38a) 


and 


Y SL ( x ,t) = A(x,t) = A 


(38b) 


Substitute the velocity profile (eq. (28a)) into equation (23a) and perform the indicated 
integrations. Then use equation (37) to obtain, for 

pMt p,n6 2 1° = - (“) + (Hi) (39a) 

0t 0X \A /j^y \A /g^ 

Fox x < X g (t) , where equations (38a) and (28b) were used in equation (23b), the follow- 
ing mass balance is obtained: 

p, — + p.n6 2 — = - (— ) (39b) 

L dt L 0X \A /j^y 


Use equation (38b) in equation (24) to obtain a mass balance for the solid layer. 



Substitute equation (33) into equation (31) to eliminate T w , and use this result along 
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with equation (38b) to evaluate the integral in equation (27) . The heat fluxes are evalu- 
ated by using equations (35), (36), and (40). This manipulation results in equation (41), 
which describes the growth of the solid layer at any position x, where x > X g (t> • 


" p s c s^ T SL ' T c> 


ty T LV ‘ T SL 


3t \A + k g b/ 


T 3A 
+ P s h SL~^ 


T SL - T c 
— + b 


(41) 


Define a subcooling parameter according to 


S 


1 “ 


c s^ T SL " T c> 


J SL 


(42) 


and rearrange equation (41) to obtain the final form of the solid layer equation. 

,2| 

o-l I / 

1 + — 1 -I 


P a L, 


s SL 


k s b 


l A + k g b j 


l 9A _ k s( T SL - T c> V T LV - T SL) 
3t A + k g b 5 


(43) 


Equation (43) relates the heat of fusion energy generated by solid layer growth (i. e. , 
freezing) and the energy that is required to cool the solid layer below the freezing tem- 
perature to the difference between the heat flux transferred to the coolant and the heat 
flux into the solid layer. 

The equation for the condensate layer, where there is no solid layer at x< X g (t) , 
is obtained next. First, substitute equation (32) into equation (30) to eliminate T w (x, t), 
and substitute that result with equations (28b) and (37) into equation (26) in order to eval- 
uate the integrals. Then, substitute equation (39b) into equation (21), and use that result 
along with equations (34) and (37) to achieve the following result: 


P, c 


ri 


(T ^- T c>^U +bk , 


i i n I 95 , n n L 5 ^ 96 P £ l /m _ m \ 9 
) + p Z L LV ^T + P Z L LV 6 — 


6" 


at 


8 


LV 


3x \ 6 + bk. 



(44) 
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Define a subcooling parameter Sg as 



c *< T LV - T c) 


J LV 


(45) 


and rearrange equation (44) to obtain the following equation, which describes the growth 
of the condensate layer where there is no solid layer at x< X (t): 

fa 


P, L 


I^IN 



bk 


1 - 


7 


l 6 + bk 



^ + Pj „L Lv a 2 ^ + ?^hfYs,-?- 


V t lv- t c> 


3x 


8 


3x \6 + bk. 


5 + bk. 


(46) 


Equation (46) relates the heat of vaporization energy liberated by condensation, the en- 
ergy required to cool the condensate below saturation, and the thermal energy convected 
away, to the heat flux energy transferred to the coolant. 

The describing equation for the condensate layer, where there is a solid layer at 
x > X g (t), is determined from equation (25). Use equations (38), (28a), and (29) to 
evaluate the integrals; then substitute equation (40) into equations (39a) and (36). Now, 
use these results with equation (38), to evaluate the heat flux terms of equation (25). 
Substitute these results into equation (25) as indicated and obtain the following equation: 


^ < T LV - T SL> f * ^ LV f + V*< T LV - T SL> f + p s L LV ~ 


at 


at 


3 m x *2 36 . „ T „«2 dd _ V T LV ' T SI^ 

+ g p l C l n(T LV " T SL )6 — + p Z L LV n6 3x ~ 5 


(47) 


Define another subcooling parameter as 


0 = c l ( t lv " T SL> 

2 L 

L LV 


(48) 


and transform equation (47) into equation (49), which describes the condensate growth 
where there is a solid layer (x > X g (t)). 
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P S L LV (S 2 ♦ l) M + Pi L LV „| S 2 + i)a 2 g 

^ Y ^ 

Coupling term 


V T LV t sl> 

6 

( 49 ) 


Equation (49) can be described the same way as equation (46), with one notable exception. 
There is a term involving 9A/3t, which is the additional heat generated by condensing 
the extra liquid that is frozen to form a layer of solid. This term couples equation (49) 
to equation (43). 

Equation (46) is solved wherever there is no solid layer (i. e. , x ^ X (t)), while 
equations (43) and (49) are solved simultaneously wherever there is a solid layer (i. e. , 
x>X g (t)). The differential equations describing the problem (i. e. , eqs. (46), (43), 
and (49)) are quasi-linear first-order simultaneous partial differential equations whose 
solution is considered in the next section. 

The solution of the problem also depends upon the point where freezing starts 
X g (t) . According to the boundary condition there (eq. (13)) the wall temperature at 
X g (t) is at the freezing point. 


T w< x=5 V l > = T SL < 50 > 

Substituting this fact into equation (32) gives an equation for the condensate thickness 6 
at x = X g <t). 


a<x s ,t> 


s 


" t sl) 


T SL - T c 


(51) 


Notice that 6_ is not a function of time or position. It is a constant for this problem 
since the parameters of equation (5 1) were assumed to be constant. In order to have 
freezing (i. e. , in order for 6 e ^ 0), the coolant temperature must be less than the 
freezing temperature (T c s Tg L ). Instead of using X g (t) as the switchover point be- 
tween equation (46) and equations (43) and (49) , it is easier to solve equation (46) for 

when 5 < 6„ and equations (43) and (49) wherever 6 > 6 . Therefore, the boundary 
s s 

conditions used in the two-phase (or condensing with no freezing) region (6 < 6 ) are 

"" s 

given by 


6<x, 0) = (P 

6(0, t) =0^ 


(52) 
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whereas in the three-phase region, where there is condensing with freezing (6 > 6 ) 

s 


6<X S , t) = 5 


s 


A<x>x s ,y = o^ 


A(X g ,t) =0 


( 53 ) 


METHODS OF SOLUTION 

The next task is to solve the partial differential equations that were derived in the 
previous section. In this section it is shown that the analytical treatment for the tran- 
sient condensing region (x < X g ) is quite different from the treatment of the three-phase 
region (x > X„). The transient condensing region can be solved by a simple exact single 
characteristic solution. Families of characteristic curves, described by ordinary dif- 
ferential equations, are necessary in the three-phase region. A direct finite difference 
solution to the partial differential equations is also employed in the three-phase region 
for comparison. The characteristic solution in the three-phase region (i.e. , families of 
characteristic curves) is somewhat unusual because one characteristic is a zero charac- 
teristic. Stability and accuracy requirements necessitated small time and position in- 
crements of either the direct finite difference or characteristic solution in the three- 
phase region, thereby resulting in long computational time. Computing time for a given 
accuracy was reduced by concentrating the smallest position elements near the start of 

freezing location X . It turned out that the stability, the accuracy, and the computation 
s 

time of the characteristic and finite difference methods were essentially equivalent. A 
single characteristic solution, which does not require lengthy computations, was pos- 
sible in the three-phase region when the term coupling the two describing differential 
equations was neglected. This term is zero only at the start of freezing and at steady 
state, so that an appreciable error might result if it is neglected. 

Details of the characteristic solutions are now discussed. This is followed by a de- 
tailed discussion of the direct finite difference solution. 
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CHARACTERISTIC EQUATIONS METHOD 
Transient Condensing Region 

As indicated in reference equation (46), the partial differential equation describing 

transient condensing for x < X c (t) , is rewritten in the following standard form. 

s 




(54) 


The result is given by 


p l L LV 


M t lv - T c> 



S 3 l bk l 


This then implies that 


2 \bkj + 6 j 


(b v^ir 


(6 + bk z ) 


n6 


S 3 n (36 + 4bk^6 


( 5 + **, i y 



l 

3x 


(55) 


dx dt 


= d6 


f 2 (6) f!<6> 

so that the single characteristic curve, described by equation (57), exists. 

.3 


nd' 


+ V ( 


36 + 4bk z )6" 


/dx\ 

v 5 > 8 ( 


w 

c f l< 6 > / S,\ So 

r l 2 


+6 


(56) 


(57) 


Assume that S 3 = 0, which is tantamount to neglecting energy storage and convection; 
then, this characteristic equation is simply given by 
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'te\ = n6 2 = " p v )g x 5 2 

dt / c 


The condensate thickness 6 is obtained from the further solution of equation (56); 
namely, 


d5 _ 1 

dx fg^) 


(59a) 


and 


d5 _ 1 

dt f 1 <5) 


(59b) 


The solutions to equations (59a) and (59b), where the boundary conditions are given by 


6(x, 0) = 0 
6(0, t) = 0 

and it is assumed that Sg = 0, are given by equations (61): 

6 = 6(x) 


or 


and 


or 




x = x(6) 


epi - 'V^lv /s 4 , bk / \ 

kA< t lv - V V 4 3 / 


6 = 6(t) 




t = t(6) 


P Z L LV 
k ^ T LV - T c> 



(60) 


(61a) 


(61b) 
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The single characteristic curve separates a steady-state solution region from a 
function-of -time -only solution region. The mathematical basis for this statement is that 
characteristic curves are curves of discontinuities in the derivatives of the dependent 
variables. Therefore, the characteristic equation (eq. (58)) determines where equa- 
tions (61) apply for S = 0. Equation (61a) describes the steady-state region, and equa- 
tion (61b) describes the function-of-time-only region, where 6 is not a function of x. 
For a given value of time t, equation (61b) is solved for 6, which is used in equa- 
tion (58) to obtain the characteristic location. Equation (61a) is then solved for 6 for 
given values of x that are less than the characteristic location. The above calculation 
procedure gives 6(x) at a given value of time t. 

The assumed constancy of the parameters in equation (51) requires that the conden- 
sate thickness at the start of freezing 6 be constant. And the single characteristic 
solution for the two-phase region means that 6 will reach 5 g uniformly with x, and 
freezing will initiate uniformly, for all values of x greater than the characteristic loca- 
tion. This characteristic location is where freezing initiates (X g (t)) in this instance. 
After the start of freezing, values of x< X (t) are in the steady-state region. There- 
fore, the point of freezing initiation X g (t) does not move after the start of freezing 

(i. e. , X (t) = X = Constant). This fixed location can be determined readily from a 
s s 

simple steady-state analysis of the condensing region. 

A further consequence of these initial freezing conditions (i.e., 6(x > X , t ) = 5 

s s s 

= Constant, and A(x > X , t ) = 0) is determined by substituting them and equation (51) 

s s 

into equation (43). From this it is learned that 0A/3t(x > X , t ) = 0. Axial conduction 

s s 

is a small effect, which is largest near X_ and tends to vanish at large x. If axial 

s 

conduction had been included in the analysis, a 9A/9x term would have appeared in 
equation (43), and 3A/3t(x > X g ,t g ) would be greater than zero near X g . 


Transient Condensing with Freezing Region 

General problem . - The partial differential equations describing the three-phase 
region x > X g (t) (eqs. (43) and (49)) can be rewritten in the following form by substitu- 
ting equation (43) into equation (49) and neglecting the subcooling terms (i. e. , s 0 and 

s 2 .0). 


36 + n5 2 - A 5 

9t 3x 6 A + Ag 


(62) 
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and 



where 



Equations (62) and (63) are put in the following matrix form: 

B^= c 

3t 3x 


where 



( 63 ) 


(64) 


(65a) 


(65b) 
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"9A" 

3x 

3 y _ 

0X 36 

_3x_ 


and 



By satisfying the requirement that the determinant 

|B - XI | = 0 


the following eigenvalues result: 


A 


1 



A 


2 



(65c) 


(65d) 


(65e) 


(65f) 


( 66 ) 


The differential equations are hyperbolic because these eigenvalues are real and unequal. 
The characteristic equations for the differential equations are 


(dt) = ° (67a) 

I / 
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I 



and 


\ = n6 2 (67b) 

n 

Because of the zero characteristic (eq. (67a)), which is a consequence of the assumption 
of no axial heat transfer, the commonly used characteristic method (e. g. , see refs. 9 
and 10) cannot be used here. Fortunately, equations (62) and (63) are recognized to be 
in the normal form of references 11 and 12. 


& + A = Constant (68) 

9t 8x 


As a consequence of this normal form, only the 8 y/ 9t derivative appears along the 
characteristic curves, and a solution can be obtained. Along the (dx/dt)^ = 0 character- 
istic, the differential equation to solve is given by equation (69) 



2 

while along the (dx/dt)jj = n8 characteristic, 

d6 = ^6 _ A 5 
dt 6 A + Ag 


(69) 


(70) 


A solution to the problem now involves the numerical solution of equations (69) and (70), 
along the characteristics given by equations (67a) and (67b), respectively. The numeri- 
cal approximations to equations (67), (69), and (70) are summarized in appendix B and 
discussed further in the section Direct Finite Difference Method. 

Coupling t erm zero. - If the coupling term 8A/ 8t of equation (49) were zero, the 
partial differential equations would be decoupled. Equation (43) would then be simply an 
ordinary differential equation to solve at each position x, and equation (49) could be 
solved by the simple single characteristic method outlined previously for the condensing 
region. The coupling term is only exactly zero at the start of freezing and at steady 
state. Following the previous single characteristic procedure, the condensate layer for 

x > X_ is described by the following equations when the coupling term is neglected, 
s 

The characteristic equation is given by 


I 
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fg(6) 

f 4 (8) 


n8 


■) 


s o 

-^+ 1 


(71) 


where 6 is obtained from the solution of equation (72) or (73). 


d6_ 1 _ ty T LV " T Sl) 


* ' 4<5> ^ L LV"{ 5 S 2 +l ) 63 


d6 


_ 1 ty T LV " T SlJ 


dt fg<6) 


^ L Lv(“ + *)< 


(72) 


(73) 


Assume that S 2 3 0; then the characteristic equation is given by 


(if)- 


n6" 


(74) 


Equations (72) and (73) are readily integrated by using the boundary conditions 


S(x, t ) = 5 = Constant 


6(X , t) = 6 = Constant 


(75) 


The result of the integration, assuming that S 2 = 0, is given by equations (76) and (77). 

x - X_ = l — 1 (76) 


s ¥ t lv- t sl>\ 4 


t-t - p ^ l lv 1 62 - a f 

s ¥ t lv- t sl> 


(77) 


Therefore, by assvuning the coupling term is zero, only algebraic equations need be 
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solved for 6 instead of partial differential equations. This would represent an appre- 
ciable computational savings. The equation for the solid layer thickness (eq. (43)) de- 
pends upon 6 at each x. Thus, it is solved as an ordinary differential equation at each 
x, where the condensate thickness is provided by equation (76) or (77). 


DIRECT FINITE DIFFERENCE METHOD 

The solution of the partial differential equations (eqs. (43) and (49)) for the 
condensing-with-freezing region by a direct finite difference method is now discussed. 
The method could also be applied to the solution of the single partial differential equation 
(eq. (46)) of the condensing region should there be different boundary conditions. 

A major problem with finite difference numerical methods is to select Ax and At 
increments such that the numerical solution is stable. This problem is considered in 
reference 13, where a method is given to calculate a time increment At, as a function 
of the Ax increment and the problem variables, which will result in a stable solution. 
Based on experience, a first-order finite difference approximation to the differential 
equations is employed. Reference 13 suggests that a backward difference be used for the 
spatial derivatives, with a forward difference for the time derivatives, and all coeffi- 
cients evaluated at the previous time. The finite difference equations used to approxi- 
mate equations (43) and (49) are found in appendix B. The following equation is obtained, 
by the method of reference 13, in order to calculate the time step required for a stable 
solution: 


At = - A - x — (78) 

2n6 2 N 

where N is an arbitrary number that is greater than 1. In order to give a constant 
value to At, so that layer profiles could be plotted at constant values of time, the terms 
of equation (78) are evaluated conservatively at the largest value of x for the problem 
(e. g. , x = L = 10 ft or 3. 05 m). This combination of finite difference approximations 
and time step resulted in a stable solution for the first time step and subsequent time 
steps. It was found that N = 1 gave a very close estimate of the maximum allowable 
At that would permit a stable solution. 

Stability does not assure accuracy. Small values of Ax and At (where At must 
always satisfy the stability criteria of eq. (78)) are chosen in order to reduce the trun- 
cation error without unduly increasing computing time (i. e. , cost) and round-off error. 
The accuracy of the solution is best determined by comparing the results of the finite 
difference calculation to results of two limiting cases, where the calculation is inher- 
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ently accurate. These two cases, which have exact solutions, are a steady-state case 
and a case where the coupling term of equation (49) is assumed to be zero. By assuming 
that the coupling term is zero, the resulting single partial differential equation for 6 
can then be solved by the simple single characteristic method with assured accuracy, 
since this is an exact solution. This analysis was worked out in the previous section. 

By also solving this same uncoupled differential equation directly by the finite difference 
method, an estimate of the finite difference method error, as compared to the equivalent 
exact solution, can be determined. In order to reduce computing time, it is desirable 
to use a small number of At and Ax increments. Figure 2 shows the percent error 
distribution along the plate, at two different times, for three different distributions of 
Ax elements, where a total of 500 Ax increments were used for each distribution. 

This figure shows that the percent error is too large in the important region of the plate 
(i.e., 0.01L<x<L) for a uniform distribution of 500 Ax increments in L = 10 feet 
(3. 05 m). The error in this region becomes acceptable if the 500 Ax increments are 
heavily concentrated near X„ (i.e., many small increments are located near X_), 
where the major changes in the phase thickness profiles occur. Distribution B of fig- 
ure 2 (150 Ax in the first 0. 05 L, 100 in the next 0. 1 L, 50 in the following 0. 1 L, and 
100 in the remainder of L) was finally used in the numerical calculations of the coupled 
case; the results of which are described in appendix B. A time increment distribution, 
based on equation (78), was used with all of these Ax distributions (i. e. , N = 1000 for the 
first 600 At increments, then N = 100 for the next 400 At, and N = 10 beyond that). 
This At distribution reduced the computation time (cost) while maintaining a stable and 


TABLE I. - COMPARISON OF STEADY-STATE RESULTS OF REFERENCE 3 AND TRANSIENT 


LAR GE - TIME a RESULTS AT x - L = 10 FEET (3. 05 m) 


Case 


Thermal resistance of plate and coolant, b, (ft 2 )(hr)(°R)/Btu; (m 2 )(K)/W 




5X10 -4 ; 0.88X10*' 4 


5X10 3 ; 0.88X10" 3 



Thickness of 

Thickness of 

Position where 

Thickness of 

Thickness of 

Position where 


solidified layer, 

condensate layer, 

freezing starts, 

solidified layer, 

condensate layer, 

freezing starts, 


A 

6 

X s 

A 

6 

X s 



Percent error - (l This r6sult ‘ Reference 3 result |\ x 1Q0 
\ Reference 3 result ) 


Uncoupled 

0. 125 

0 . 1 

3 

0 . 1 

0 . 1 

0.6 

transient, S = 0 







Coupled 

.6 

.6 

3 

.5 

. 1 

.6 

transient, S = 0 








a In the transient cases, steady state was considered to be when no change occurred in the third significant figures of 
A and 6 at x = 10 ft (3. 05 m) for 500 calculations. 
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accurate solution (N > 1) . Table I contains a comparison of the approximate steady- 
state results from the finite difference calculation (i. e. , d/dt — 0) and exact results 
from the closed-form steady-state analysis of reference 3. The agreement is good. 

The characteristic solution of the coupled partial differential equations for the 
three-phase region was necessarily also solved by an approximate finite difference 
method. The accuracy and stability of this numerical solution are now compared to the 
previous direct finite difference solution results. One comparable case was computed, 
where the numerical approximations were of the same order and the parameters, Ax 
and At increments, and stability criteria were the same. The computed phase thick- 
nesses were within 1/2 percent, and the two methods required very nearly the same com 
puting time to reach the same problem time. It was found that equation (78), with N = 1 
also gave a close estimate of the maximum allowable At for a stable solution for the 
coupled characteristic equations. 


RESULTS AND DISCUSSION 


In this section the results of the transient analyses are discussed. The first part of 
this section is devoted to transient condensing of a pure vapor on a nonisothermal in- 
clined plate where there is no freezing. The second part takes up transient condensing 
of a pure vapor with freezing. 


TRANSIENT CONDENSING OF A PURE VAPOR ON 
A NONISOTHERMAL INCLINED PLATE 

The three-phase problem is made up of a two-phase region and a three-phase re- 
gion. The problem of transient condensing on a nonisothermal inclined plate is the two- 

phase part of the whole solution. If the plate is too short to have freezing (i. e. X > L) 

s 

or T c >T SL’ only the two-phase part of the whole problem need be considered. If it is 
further assumed that b = 0, then the special case of transient condensing on an isother- 
mal plate results. This case was considered in reference 5. From the steady-state re- 
sults of reference 3, which are discussed in appendix C, it is known that freezing is dif- 
ficult to achieve during the condensation of a pure vapor on a vertical plate, no matter 
how cold the coolant. Generally, very low values of pressure or coolant thermal resist- 
ance, or a long plate inclined more toward the horizontal, or some small concentration 
of a noncondensable gas is necessary. Therefore, it can be expected that transient con- 
densing with no freezing should be a more common occurrence for a vertical plate than 
condensing with freezing. 
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A single characteristic curve describes the two-phase condensing region. This 
curve, which is described by equation (57) (or, when S = 0, by eq. (58)), separates two so- 
lution regions in the x, t plane. In one region steady-state condensing occurs, and in the 
other region condensing is only a function of time. Figure 3 shows a group of such single 
characteristic curves for water condensed on a vertical plate at a given pressure and 
coolant temperature. The coolant and plate thermal resistance b is varied as a param- 
eter. The characteristic curve indicates the time required to attain steady-state con- 
densing at various positions x for a particular set of conditions. For the case consid- 
ered in figure 3, a good representative time to attain steady-state for a 1-foot- (0. 3-m-) 
long plate would be about 1 second. From equations (58) and (61b) or figure 3, it can be 
seen that it takes longer to attain steady state, at a given x, as the coolant and plate 
thermal resistance b increases or as the plate becomes more horizontal (g — 0) . The 
effect of an increase in coolant temperature or pressure results in a much smaller in- 
crease in the time to attain steady state. The broken lines in figure 3 were determined 
from equations (57) and (59b) for the case where Sg ± 0 in order to determine its effect. 
Clearly, the effect of subcooling could be neglected (i. e. , Sg = 0). 

TRANSIENT CONDENSING WITH FREEZING ON 
A N0NIS0THERMAL INCLINED PLATE 

In this section, results of the less common but interesting situation of transient 
condensing with freezing are considered. This three-phase situation would tend to 
occur in practice when T c < T SL’ and the pressure is low, the plate is nearly hori- 
zontal, or a small amount of noncondensable gas is present in the vapor. The additional 
complication caused by the inclusion of the noncondensable gas effect in a transient situ- 
ation is beyond the scope of this report; however, it is analyzed for steady state in ap- 
pendix C. 

It was previously shown that a simple single characteristic solution would result only 
when the coupling term in equation (49) was neglected. This term, which accounts for 
the heat of vaporization liberated by the additional condensate mass needed to form the 
growing solid layer, is exactly zero only at the start of freezing and at steady state. 
During the transient process, the numerical value of the coupling term is about the same 
size as the 96/8t term in equation (49); therefore, it cannot generally be neglected in 
spite of the desirable simplification that results. The error caused by neglecting the 
coupling term is clearly shown by comparing the coupled (numerical solution) and uncou- 
pled (single characteristic solution where there is no coupling term included) results of 
figure 4. Figures 4(a) and (b) display the growth rate of the condensate and solid layers 
using the coupled and uncoupled analyses at two positions on the plate (x = 5 ft (1. 5 m) 
and 10 ft (3m) for b = 5xl0 -3 and 5X10" 4 (ftfyhr^Rj/Btu (0. 88X10 -3 and 0. 88X10" 4 
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(m^Kj/W), respectively). The coupling term, as shown in figure 4, retards the attain- 
ment of steady state. In the section, METHODS OF SOLUTION it was shown that 
3A/3t (x > X c ,t ) = 0 because axial conduction was neglected. The zero derivative is 
lost on figures 4(a) and (b) because of the scale of the plot. Figure 5(a) contains the 
phase thickness profiles (S(x), A(x)) at various times for the coupled solution where 
b = 5xl0 -3 (ft 2 ) (hr) (°R) /Btu (0. 88xl0 -3 (m^Kj/W). Figure 5(b) contains a similar com- 
parison for a lower value of b. 

For a chronological description of the phase growth consider figure 5(a), which 
shows the phase growth history for a coupled solution where b = 5x10" 3 (ft 2 ) (hr) (R) /Btu 
(0. 88xl0” 3 (m^K^/W) . Initially (t = 0), the plate is cold but dry (6 = 0). The conden- 
sate layer starts to grow an instant later and with increasing time tends to ’’fill up” the 

steady-state profile until t = t when 6 reaches the condensate thickness required to 

s 

start freezing 6 . From this point in time on, there exists a two-phase region above 
s 

the fixed-start-of -freezing location at X g , and a three-phase region below that point. 

The two-phase region, having reached steady state when freezing started, does not grow 
anymore. In the three-phase region, the condensate and solid layers continue to grow 
until they fill up their steady-state thickness profiles. Essentially, the same description 
applies to figure 5(b), where b = 5X10" 4 (ft 2 ) (hr) (°R) /Btu (0. 88xl0" 4 (m 2 )(K)/W). How- 
ever, freezing starts very near x = 0 and the phases grow thicker. 

Consider the special case of b = 0, which results in a constant wall temperature 
(T = T c = Constant). According to equation (51), 6 g = 0, so that freezing starts imme- 
diately at t = 0 and all along the plate. This academic case cannot be solved by the nu- 
merical methods discussed but can be very closely approximated by a case where 
b = 10”® (ft 2 ) (hr) (°R) /Btu (0. 18x10”® (m^Kj/W. The computing time is unfortunately 
very long for such a low value of b. 

In spite of the inaccuracy of the uncoupled solution, this simple solution can give a 
crude estimate of the time to attain steady state. The error increases as b decreases. 
For example, according to figure 6 steady state is attained on a 1-foot (0. 3-m) vertical 
plate in about 1 second. 


CONCLUDING REMARKS 

Transient condensing on an initially dry inclined plate of no thermal capacity that 
is cooled by a fluid below the freezing point of the vapor was considered in this anal- 
ysis. The problem was characterized by a two-phase condensing region and a three- 
phase region where the condensate freezes to form a growing solid layer. A Karman- 
Pohlhausen integral method, with linear temperature profiles across the phases was em- 
ployed. This reduced the problem to the solution of simultaneous quasi-linear partial 
differential equations, where the condensate and solid thicknesses were the dependent 
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variables. A single characteristic curve, represented by algebraic equations, de- 
scribed the transient condensing (two phase) region. Families of characteristic curves, 
represented by ordinary differential equations, were involved in the solution for the 
three-phase region. A solution involving a finite difference equation approximation to 
the partial differential equations was also employed in the three-phase region for com- 
parison. Small time and position increments were required for the numerical solution 
by either the finite difference or the characteristic method in the three-phase region in 
order to achieve an accurate and stable computation. This resulted in a long computing 
time. The stability, accuracy, and computation time of the characteristic and finite dif- 
ference methods were found to be essentially equivalent. A single characteristic solu- 
tion, which required no long computation, was possible in the three-phase region only 
when the term coupling the two describing partial differential equations was neglected. 
This coupling term accounts for the heat liberated by the extra liquid that must be con- 
densed in order to form the growing solid layer. An appreciable error can occur when 
this coupling term is neglected. This term is negligible only at the start of freezing and 
at steady state. 

The results of this analysis indicated that freezing starts at and beyond a particular 
point on the cold plate as soon as the insulating condensate layer is thick enough for the 
wall temperature to drop to the freezing temperature at that point. The location of the 
start of freezing does not change after freezing begins. This fixed location is a conse- 
quence of the single characteristic solution for the transient condensing region. The 
solid and condensate layers, which have a shape similar to boundary layers, were found 
to grow until they attain a steady-state profile. 

The analysis of transient condensing with freezing was shown to specialize to a 
number of interesting cases which are discussed in appendix C. Some examples are 
transient condensing and freezing on a horizontal plate, transient condensing on a non- 
isothermal inclined plate, and steady-state condensing and freezing (or melting) on a 
nonisothermal inclined plate. With minor changes in the initial conditions, the analysis 
could handle transient condensing with melting; and with a minor modification, transient 
freezing of a liquid flowing past a cold plate could be adequately analyzed. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 26, 1969, 

120-27-04-27-22. 
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APPENDIX A 


SYMBOLS 


A 

b 


V 



c 



*!<«>, *2<°> 


Sc 

So 

g x 


h 

I 

k 

L 

^LV’ L LV < ' P ^ 

l sl 

(m/A) c d> (m/A) fr 
(m/A) LV 


(m/A) 


SL 


N 


n 


P 

(Q/A) c 

q 

s l’ s 2’ s 3 


constants, defined by eq. (64) 

thermal resistance of plate and coolant, b = 1/h + d A , 

(ft; (hr) (°R) /Btu; (nO(K)/W 

specific heat at constant pressure, Btu/(lbm)(°F); J/(kg)(K) 
thickness of plate, ft; m 
functions of 6 

o p 

standard conversion factor, ft/hr ; m/sea 

9 9 

steady acceleration of gravity, ft/hr ; m/sec^ 

steady body force acceleration along inclined plate, g^= g Q sin <p, 
ft/hr m/sec ^ 

surface heat-transfer coefficient, Btu/(ft^)(hr)(°R); W/(m^)(K) 
unit diagonal matrix 

thermal conductivity, Btu/(ft)(hr)(°R); W/(m)(K) 

length of plate, ft; m 

heat of vaporization, Btu/lbm; J/kg 

heat of fusion, Btu/lbm; J/kg 

mass flux (condensing and freezing), lbm/(hr)(ft^); kg/(m^)(sec) 

mass flux across liquid-vapor interface in y-direction, 
lbm/(hr)(ft 2 ); kg/(m^)(sec) 

mass flux across liquid-solid interface in y-direction, 
lbm/(hr)(ft 2 ); kgAmfysec) 

number greater than 1 

[(fy - P v )S x ]/m 1 , 1/ (ft) (hr) ; l/(m)(sec) 

pressure, psia; N/m^ 

heat flux to coolant, Btu/(ft^)(hr); W/m^ 

heat flux, positive in y-direction, Btu/(ft^)(hr); W/m^ 

subcooling parameters, defined by eqs. (42), (48), and (45) 
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T i 


LV 

r SL 




At. 


U, 


u. 


W 


x 

X A’ X R’ X C 

Ax i 

(dx/dt) p (dx/dt) c 

Y SL’ Y SL< x > ^ 
Y LV’ Y LV^ X ’ 

y 

Y 

A, A(x), A<x, t) 

5, 6(x), 6(x, t> 

6 s 

M 

P 

<P 


temperature, °R; K 
interface temperature, °R; K 
boiling temperature, °R; K 
freezing temperature, °R; K 

superheated vapor temperature far from plate, °R; K 

wall temperature (temperature of the surface of plate in contact 
with phases), °R; K 

time from start of condensation, hr; sec 

nonconstant time increment, defined in appendix B, hr; sec 

time when freezing starts, hr; sec 

bulk vapor velocity, ft/hr; m/sec 

velocity of liquid in x-direction relative to fixed axes, ft/hr; m/sec 

velocity of liquid in y-direction, relative to fixed axes, ft/hr; 
m/sec 

noncondensable gas mass fraction 

position along plate where freezing starts, ft; m 

x-coordinate direction along plate from start of cold section, ft; m 

points in mesh, defined in appendix B 

nonconstant x increment, defined in appendix B, ft; m 

slope of characteristic curves 

location of melt surface (SL) from plate, ft; m 

location of condensing surface (LV), ft; m 

y-coordinate, normal to plate, ft; m 

term defined by eq. (65) 

thickness of solidified layer, ft; m 

thickness of condensate layer, ft; m 

thickness of condensate layer where freezing starts, ft; m 

dynamic viscosity, lbm/(ft) (hr) ; N-sec/m 2 

3 3 

mass density, lbm/ft ; kg/m 
angle of plate from horizontal, rad 
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Subscripts: 


c coolant 

cd condensing 

fr freezing 

i indicates x-direction in numerical calculations 

j indicates t-direction in numerical calculations 

l liquid 

m plate 

s solid 

v vapor 
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APPENDIX B 

DISCUSSION OF NUMERICAL CALCULATIONS 


The numerical calculations for the finite difference and characteristic solutions of 
the three-phase region are discussed further in this appendix. In particular, the finite 
difference approximations to the coupled partial differential equations in the text are de- 
scribed. 


CHARACTERISTIC CALCULATION 


The characteristic solution of the three-phase region involves the solution of equa- 
tion (69) along characteristic I, given by equation (67a), while equation (70) is solved 
along characteristic n, given by equation (67b) . Consider figure 7, which shows typical 
characteric curves for this problem and the grid used for the calculation. The values of 
6. . (6 at point B), 6 at point R, and x R are to be determined. To solve equation (70) 
along the II characteristic, given by equation (67b), linear approximations to these equa- 
tions between points B and R are used. The value of x at point R is determined from 
the following approximation to equation (67b): 


X R “ X C 



(79) 


while 6. • 
J 


is determined by the following approximation to equation (70) 




A R + A 2/ 


(80) 


Values of 6^ and Aj^ are determined by the following linear interpolations at x^ be- 


tween values of 6 and A at x^ and x^. 




At- . 


At- 




(81a) 
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I 


and 


A R=< A i,j-l> 


At. 

l-n-±(6. . ,V 
AXj V 1,3-1/ 


At. 


+ n — i- (s. . * r(A. 1 . 1 ) 
Av \ 1,3-1/ v i-l, J-l 7 


AXj 


(81b) 


Since the intervals are not constant, AL = t^ - t. ^ and Ajc = - x i Now the value 

of A at point B A^ j is determined by solving equation (69) along the I characteristic. 
This solution simply involves the solution of a quadratic equation in Aj for the largest 
positive root, making use of the previously calculated values of 6. .. ’ 


i,3 



A 2 A i,j-l+ M j A 4fl- / 'j 

\ i, 3 / 

L \ i, 3 /_ 


= 0 (82) 


FINITE DIFFERENCE CALCULATION 

The direct finite difference solution to equations (62) and (63) involves an explicit 
solution where a first-order forward difference is used for the time derivatives and a 
first-order backward difference is used for the x derivatives (see fig. 8). 


6. . « - 6. 


at 

7-7 7 w 

Al J 

3A _ 

A i,j+1 " A i, j 

at 

At i 

36^ 

6. . - 6. 

1,3 1-1,3 

ax 

Ax- 


l 


(83) 


J 


The coefficients in the partial differential equations that involve the dependent variables 
are evaluated as 6 and A. .. Equation (62) is approximated by the following differ- 


ence equation which is solved explicitly for the unknown 6. . 1 . 


6. . 1 = 6. . + At. 
i,J+l i,3 3 


Ki A i,3 + A 2 >J \ ^i )\ 


(84) 
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Equation (63) is similarly approximated, resulting in a quadratic equation for A i ,. +1 , 
where 6. ^ is supplied from equation (84) and the largest positive root is used! 


/ \2 / At i A 

°-( A i,j+l) + A i,j+1^2 + -^— 


A 2 A i, j + At j A 4( A l " 


A 2 A 3 



(85) 


This numerical calculation procedure will not start up at the start of freezing (i. e. , 
at t = t g ). Fortunately, the coupling term is zero at this time so that the single char- 
acteristic approximation can be used in calculating the first time step. Once the first 
step is calculated, the numerical procedure is used. 
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APPENDIX C 


DISCUSSION OF SPEC IAL CASES 

An interesting feature of the three-phase problem is the wealth of special cases that 
can result by simplifications of, or minor modifications to, the partial differential equa- 
tions. Some of the special cases have not appeared in the literature. 

Two of the special cases have been discussed in the main text; they are the three- 
phase transient problem, where the inclined plate is isothermal, and transient condens- 
ing on a nonisothermal inclined plate. Two special cases that have been analyzed in 
reference 3 are summarized here because the results give the reader a more complete 
understanding of the physics of the problem. These cases are transient condensing and 
freezing on a horizontal plate and steady-state condensing and freezing (or melting) of a 
pure vapor on an inclined plate. The important effect of a noncondensable gas on the 
steady-state problem was also considered in reference 3; it too is summarized in this 
section. 

Three other transient special cases are discussed that require minor modifications 
to the describing partial differential equations: melting of a thin slab of solid by its 
vapor, freezing of a liquid that flows past a cold plate, and finally the three-phase prob- 
lem when the coolant and vapor parameters (e. g. , T , h , and P) are known functions 
of time and/or position. No numerical results are worked out for these three cases. A 
major modification is required for a case where there is surface evaporation, or where 
the plate has finite thermal capacity. 


STEADY-STATE CONDENSING WITH FREEZING 
(MELTING) ON A NONISOTHERMAL PLATE 

In this section, the special case of steady state is considered, and the effect of 
parameters upon condensing with freezing is demonstrated in detail. This special case 
has been considered in reference 3 and is summarized herein in order to give a more 
complete physical discussion of the three-phase problem. By neglecting the time deriv- 
ative terms and energy convection (i. e. , Sg s 0, Sg s 0), equations (46), (43), and (49) 
directly simplify to the steady-state equations of reference 3, which describe steady- 
state condensing of a pure vapor on a nonisothermal inclined plate. These equations are 
manipulated and solved to show the conditions for freezing. The effect of noncondensable 
gas on freezing is also demonstrated by a minor modification to the steady-state anal- 
ysis. 
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Pure Vapor 


Setting the time derivative terms equal to zero, neglecting energy convection (i. e. , 
Sg = 0, Sg = 0), and integrating according to the boundary conditions (eqs. (52) and (53)) 
results in the following steady-state equations for the condensate liquid and solid layer 
profiles. Where there is no solid layer, at x< X . 

- T C>* (86) 

4k l 3 - OyK 

Where there is a solid layer, at x > X , the following are the condensate and solid layer 

fa 

thickness profile equations: 



m ^^ T LV ~ T SL^ x ~ X s> 
p Z L LV^Z “ p v )g x 


and 


A = 



(87) 


( 88 ) 


where 6 g is given by equation (51) and the wall temperatures by equations (32) and (33). 

According to equation (51), these results simplify further for a case where the con- 
densate is warmer than the freezing point (i. e. , T c > Tg^), since no freezing is possible 
(6 s < 0). In this case, the remaining equation (eq. (86)) describes steady-state condens- 
ing on a nonisothermal inclined plate. If it is further assumed that b s 0, so that the 
wall temperature is constant, the equation is exactly the one solved by Nusselt (ref. 6) . 

Equations (32), (33), (51), and (86) to (88) were solved for T , 6 , 6 , and A to 

w s s 

indicate the effect of the system parameters for a water system. Figure 8(a) contains 
the nominal case from which one parameter is varied for each subsequent part of fig- 
ure 8. From figures 8(a) to (d), it is clear that an increase in coolant temperature, 
pressure, and/or coolant thermal resistance b reduces the condensate thickness by a 
small amount. However, this small change can greatly reduce the solid thickness and 
the possibility of having a solid form on a given plate. 

The heat flux in all cases of figure 8 is of the order of 10 & Btu/(ft 2 )(hr)(3. 15X10 5 
W/m 2 ), which is above the "burnout flux” for many cryogenic coolants (e. g. , see 
ref. 14). Therefore, the cases of figure 8 would be limited to a value of b correspond- 
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ing to film boiling (e. g. , b < 0. 01 (ft^hr^Rj/Btu or 1. 8xl0“ 3 (m^lp/W) in many real 
cases. For an example, consider a liquid-nitrogen-cooled vertical plate at the condi- 
tions of figure 8(d); clearly, there is no freezing on the plate. A considerable reduction 
in pressure, or a reorientation of the plate toward the horizontal would again result in 
freezing on the plate for the conditions of figure 8(d) . Therefore, in summary, the cool- 
ant thermal resistance b would not, generally, be small enough in practice to allow a 
solid layer to form on a short vertical plate from the pure slowly moving vapor, no mat- 
ter how cold the coolant, unless the pressure is very low or the plate is nearly horizontal. 

Effect of a Small Quantity of Noncondensable Gas 

If a small quantity of noncondensable gas (e. g. , air) is present in the vapor (e. g. , 
steam), the air would concentrate at the liquid-vapor interface because of the slow dif- 
fusion of air back into the bulk vapor. The air concentration would lower the partial 
pressure of the water vapor (P^ < P) at the liquid-vapor interface, thereby lowering the 
liquid-vapor interface temperature T^(P.) below saturation T^y This decreases the 
driving temperature difference across the condensate layer T. - T w and effectively in- 
creases the insulating ability of the condensate layer. Reference 15 has analyzed the 
effect of a noncondensable gas (air) on condensing of a vapor (steam) upon a constant 
temperature vertical cold plate (T^ = Constant). The results of the analysis for air in 
steam can be adapted to this analysis to determine Tj. The wall of reference 15 is es- 
sentially replaced here by the solidified layer where T' = T cx = constant. With T. 
known, as a function of the total pressure P and the air mass fraction W, the effect of 
W upon 6 g , X g , 6, and A can be determined by replacing T LV by T. in equations 
(51) and (86) to (88), respectively. By using the curves from reference 15, values of 
T. are determined for W = 0 (pure vapor), 0.01, and 0. 05 and P = 1. 7 and 14. 7 psia 


TABLE n. - INTERFACE TEMPERATURE 
CALCULATED FROM REFERENCE 15 


Total 

pressure, 

P 

Noncondensable gas mass fraction, W 

(Ti=T LV ) 

0.01 

0.05 

Interface temperature, T. 

psia 

N/m 2 

°F 

K 

°F 

K 

°F 

K 

14. 7 
1.7 

10 5 

1. 17X10 4 

212 

120 

373 

322 

92. 5 
44.0 

307 

280 

65. 5 
38.5 

292 

278 
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(11. 7x10^ and 10^ N/m^). These values are listed in table II. With these values of T^, 
the phase thicknesses shown in figure 9 are obtained. Clearly, a small fraction of non- 
condensable gas greatly increases the solid thickness and the chance of having a solid 
layer at a given position x' . 


TRANSIENT CONDENSING WITH FREEZING ON A HORIZONTAL COLD PLATE 

Far out along the inclined plate, or for a horizontal plate, where the 9/9x terms 
are negligible, the problem becomes a function of time only. The equations describing 
this case are obtained by setting 9/9x b 0 in the governing partial differential equations 
(eqs. (43), (46), and (49)). These resulting equations are the same as those found in 
reference 3, where the Karman-Pohlhausen method was used. Figure 10 contains a plot 
from reference 3 of the phase growth rates for various conditions. Figure 10 shows 
that for b > 0 a condensate layer grows until it is thick enough for the wall temperature 
to reach the freezing point, at which time a solid layer starts to grow. As b approaches 
zero, this time lag goes to zero such that at b = 0 the solid and condensate both start 


TABLE HI. - ERROR IN CALCULATION OF CONDENSING SURFACE LOCATION BY 
VARIOUS METHODS FOR ONE -DIMENSIONAL CASE 


WHERE T w = CONSTANT (OR b = 0) 


Time, 

min 


0. 06 
.3 
.6 

3.0 

6.0 
30. 0 
60.0 


Exact solution for 

superheated 

a 

vapor 


Karman-Pohlhausen method with linear temperature profile, 
saturated vapor, and nearly zero thermal resistance,^ 


S * 0 


All S = 0 


Percent error in Y L y compared to exact solution 
for a saturated vapor 


2. 5 


1. 5 
1. 3 


1 . 2 
1 . 2 


a Superheat = 125° F (70 K). 

Thermal resistance, b = 10~ 6 (ft 2 )(hr)(°R)/Btu or 0. 18xl0 -6 (m 2 )(K)/W is an 

adequate approximation to b = o, and is computationally necessary. The 

wall temperature T is constant for b = 0. This limiting case is 

equivalent to the exact solution where T = Constant. 

w 
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growing at the same time. The growth rates of both phases, for b = 0 and for large 
time at any value of b, are proportional to (time) The results of these one- 
dimensional equations were compared in table in (taken from ref. 3) for the special case 
of b = 0(i.e., a nearly constant wall temperature) to the corresponding constant wall 
temperature exact solution. This comparison shows that the approximate Karman- 
Pohlhausen integral method with S * 0 gave results that are very close to the exact so- 
lution results. The comparison in table III also indicates that even a greatly superheated 
vapor has little effect. 


NEGLIGIBLE SUBCOOLING 

By neglecting the energy storage and thermal energy convection terms in the energy 
equation in comparison to the conduction terms, a considerable analytical simplification 
is possible. By this assumption, equation (3), for example, reduces to 

9 2 T, 

= 0 (89) 

3y 2 

which results in a linear temperature profile. This type of assumption is used for each 
phase layer and is generally restricted to cases where the phase layers are thin and 
there is a phase change. Carrying out the analysis as before, with equation (89) replac- 
ing equation (3), results in differential equations for the condensate layer growth that 
are much more easily derived. The same equations can also be derived from equations 
(43), (46), and (49) by setting Sj, Sg, and Sg equal to zero. 

Neglecting all the subcooling terms (S s 0) can produce a significant error, as shown 
by table III; however, if the subcooling in the solid layer (Sj) is not neglected, this er- 
ror becomes acceptable. 


MELTING OF A THIN SLAB OF SOLID BY ITS VAPOR 

Consider a thin slab of solid (e. g. , ice) of initial thickness A(x, 0) that is attached 
to an inclined plate where the coolant is at constant T c and b. Saturated steam is 
suddenly introduced, and the ice melts while the vapor condenses and the liquid rims 
down over the plate and ice layer. The ice melts most rapidly at the leading edge where 
the insulating condensate layer is the thinnest. In time, steady state will be attained, 
where there will be a condensate layer and perhaps a solidified layer over some portion 
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of the plate. For the present analysis to fit this physical case, the initial solid layer 
must be thin so that the time lag required to heat the solid to a linear temperature pro- 
file is small compared to the time required to reach steady state. The previous analysis 
(eqs. (43), (46), and (49)) can inherently handle melting as well as freezing. Whether 
there is melting or freezing depends upon the boundary and initial conditions compared 
to steady state. Equation (43) indicates that melting (i. e. , dA/dt is negative) will occur 
if the initial ice thickness is greater than that of steady state for those same conditions. 

If the initial thickness is smaller than that demanded by steady state, freezing will 
occur. The same differential equations and numerical program will therefore handle 
changes from one steady state to another, requiring only minor modifications, such as in 
the initial conditions. 


TRANSIENT FREEZING OF A FLOWING LIQUID UPON A COLD PLATE 


Consider a liquid flowing over a cold flat plate where conditions suddenly change 

somewhat such that a solid layer begins to form and grow on the plate from the liquid. 

Equation (43) for the solid layer, which was derived on the assumption that T (y) is 

s 

nearly linear, is adapted to this problem by replacing the condensate layer thermal re- 
sistance 6/k by an equivalent resistance for the flowing liquid l/h^; and T^y is re- 
placed by T^ as indicated by equation (90) . 


*s W “ yl 


k s b 




lk g b + A, 


- 1 


^=-V t i- t sl> + 


T SL - T c 

A + b 


(90) 


There are considerable difficulties in estimating for this transient turbulent liquid 
flow, where mass is removed from the boundary layer (i. e. , suction) to form the grow- 
ing solid layer. It is therefore probable that the approximations concerned with the 
solid layer that were needed to obtain equation (43) are of relatively minor significance. 
The simplest assumption is that h^ is time -invariant. Indeed, the experimental study 
of reference 1 indicated that h^ was approximately time-invariant, at least at large x. 
Theoretical or empirical results for flow past a plate could then be used to estimate h^. 
Since h^ is a function of position, decreasing with x, the solid layer will have the 
shape of a boundary layer. The solution of this problem involves solving the above ordi- 
nary differential equation (eq. (90)) in time at a number of x positions, beyond where 
freezing starts at X . The parameters T., T , and b are described by known con- 
stants, differential equations, or functions. The start of freezing X g is determined 
from the location where the wall temperature reaches the freezing point. 
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V- X s>= T SL= . < 91 > 

hli l 

NONCONSTANT PARAMETERS 

In the formation of the three-phase analytical model, it was assumed that the coolant 
thermal resistance b, coolant temperature T , and system pressure P were held con- 
stant during the transient process. It was also assumed that the plate had negligible 
thermal capacity. These assumptions resulted in a single direction of mass flow (e. g. , 
freezing, but not freezing followed by melting). By relaxing the constant parameter re- 
strictions such that the parameters P, T , and b are known functions of time and posi- 
tion, it is possible to have a nonconstant X , and also freezing at one time and place 
and melting elsewhere. There is no inherent reason why the numerical program for the 
differential equations could not be simply modified to handle this case of known functions 
of P(x, t), T c (x, t), and b(x, t). Accounting for the thermal storage of the plate, even 
where it is assumed that the plate is a simple thermal capacitor of no thermal resist- 
ance, would require major changes in the present analysis. 
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Figure 1. - Stages of phase layer growth on vertical cold plate where coolant temperature 
is below freezing temperature (T c <T SL ). 
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Figure 2. - Changes in percent error distribution caused by changing concentration 
of total of 500 small and large Ax increments over 10-foot- (3. 05-m-) long plate. 
Position where freezing starts, X $ = 3.82x10"^ feet (1. 16xl0“ 3 m). 
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Time from start of condensation, t, sec 



io " 5 io ' 4 icr 3 io " 2 10' 1 io° 10 1 


Position on plate, x, m 

Figure 3. - Single characteristic curves for transient condensing of water on a 
nonisothermal vertical plate. Pressure, 14. 7 psia (IO 5 N/m 2 ); coolant tem- 
perature, 140° R (77 K). 
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Condensate (6) and solid layer (A) thicknesses, 


24 X io 4 Position on plate. 



Time from start of condensation, t, hr 


0 10 20 30 

Time from start of condensation, t, sec 

(a) Thermal resistance of plate and coolant, 5xl0“ 3 (ft 2 Khr){°R)/Btu (0. 88xl0“ 3 (m 2 )(K)/W). 



Time from start of condensation, t, hr 


0 10 20 30 40 50 60 

Time from start of condensation, t, sec 

(b) Same as 4(a) but thermal resistance of plate and coolant, 5x10 " 4 (ft 2 )(hr)(°R)/Btu 
(0. 88x10 " 4 (m2)(K)/W). 

Figure 4. - Comparison of phase growth histories at 5 feet (1.53 m) and 10 feet (3. 05 m) along plate, 
by coupled and uncoupled computations. Water; ail S 0 0; pressure, 14.7 psia (10 5 N/m2) ; 
coolant temperature, 140° R (77 K). 


47 



Time from start of 



Position on plate, x, m 


(a) Thermal resistance, b= 5xl0~ 3 (ft 2 )(hr)(°R)/Btu (0. 88x10 ~ 3 (m 2 )(K)/W). 

Time from start of 
condensation, 
t, 



0 12 3 

Position on plate, x, m 


(b) Same as 5(a) but thermal resistance, b = 5x10 4 (ft 2 )(hr)(°R)/(Btu 
(0. 88x10-4 {m2){K)/W). 

Figure 5. - Coupled calculation of condensate and solid layer thickness profiles at various times 
after start of condensing of vapor upon initially dry vertical cold plate. All subcooling param- 
eters S = 0; pressure, 14. 7 psia (10 5 N/m 2 ); coolant temperature, 140° R (77 K). 
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Figure 6. - Uncoupled approximation in three-phase region (x> X s ) resulting in single characteristic 
curves. Subcooling parameter $2 = 0; pressure, 14. 7 psia (HP 5 N/m^); coolant temperature, 140° R 
(77 K). 
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(a) Characteristic calculation grid. 
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(b) Direct finite difference calculation grid. 


Figure 7. - Numerical calculation grids. 
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(a) Nominal case: pressure, 14.7 psia (10 5 N/m 2 ) ; 
coolant temperature (liquid-nitrogen coolant), 
140° R (77 K); thermal resistance of plate and 
coolant (nucleate boiling coolant), 5 xl 0 ~ 4 
(ft 2 )(hr)(°R)/Btu (0. 88xl0'4 ( m 2)(K)/W); bulk 
vapor velocity, 0; subcooling parameters, S c 0. 




(c) Nominal case, except pressure is 50 psia 
(3.45xl0 5 N/m 2 ). 
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(d) Nominal case, except thermal resistance of plate 
coolant (film boiling coolant) is 0.01 (ft 2 )(hr)(°R)/ 
Btu (0. 18x10-2 (m2)(K)/W). 
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Figure 8. - Steady-state condensing and freezing 
of pure water vapor on cold vertical plate for a 
nominal case and variations. 
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Figure 10. - Phase layer growth history as water vapor 
condenses and freezes on initially dry horizontal cold 
plate. Pressure, 14. 7 psia (10 5 N/m*) ; coolant temper- 
ature, 140° R (77 K). 
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(a) Same case as figure 8(d), W * 0 (pure vapor). 



Position on plate, x, m 
(b) W-0.05. 


Figure 9. - Effect of noncondensable gas (air) on steady-state 
condensing and freezing. Same conditions as in figure 8(d) 
but for various noncondensable gas mass fractions W. 
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